%% Set paths:

clear;clc;close all

%% Read data

data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv');
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

%% Table OA.6, Panel A) $10 Billion Threshold (Pre and Post Dodd-Frank) (avg of past 4 quarters assets):

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets_past_4_qtr_average >= amin & data_orig.assets_past_4_qtr_average <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; %data(data.post == 0,:).assets_past_4_qtr_average;
apost = data(data.post == 2,:).assets_past_4_qtr_average;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets_past_4_qtr_average')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets_past_4_qtr_average')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_avg4qtr.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');


%% Table OA.6, Panel B) $10 Billion Threshold (Pre and Post Dodd-Frank) (past calendar):

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets_past_calendar_year >= amin & data_orig.assets_past_calendar_year <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; %data(data.post == 0,:).assets_past_calendar_year;
apost = data(data.post == 2,:).assets_past_calendar_year;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets_past_calendar_year')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets_past_calendar_year')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_pastCldDate.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');



%% Table OA.6, Panel C) $10 Billion Threshold (Pre and Post Dodd-Frank) (past min value):

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets_past_4_qtr_min >= amin & data_orig.assets_past_4_qtr_min <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; %data(data.post == 0,:).assets_past_4_qtr_min;
apost = data(data.post == 2,:).assets_past_4_qtr_min;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets_past_4_qtr_min')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets_past_4_qtr_min')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_pastMinValue.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');




